MODULE HVAC_ROUTINES
 
! Compute the HVAC mass and energy transport
 
USE PRECISION_PARAMETERS
USE GLOBAL_CONSTANTS
USE MESH_POINTERS
USE DEVICE_VARIABLES
USE CONTROL_VARIABLES
USE COMP_FUNCTIONS, ONLY: SECOND, CHECKREAD, SHUTDOWN
USE MEMORY_FUNCTIONS, ONLY: ChkMemErr
 
IMPLICIT NONE
PRIVATE
CHARACTER(255), PARAMETER :: hvacid='$Id: wall.f90 4624 2009-08-28 13:02:16Z randy.mcdermott $'
CHARACTER(255), PARAMETER :: hvacrev='$Revision: 4624 $'
CHARACTER(255), PARAMETER :: hvacdate='$Date: 2009-08-28 09:02:16 -0400 (Fri, 28 Aug 2009) $'
CHARACTER(100) :: MESSAGE
REAL(EB), ALLOCATABLE, DIMENSION(:,:) :: LHS
REAL(EB), ALLOCATABLE, DIMENSION(:) :: RHS
INTEGER :: ITER,ITER_MAX=10
PUBLIC HVAC_CALC,GET_REV_hvac,READ_HVAC,HVAC_BC_IN,FIND_NETWORKS
 
CONTAINS

SUBROUTINE READ_HVAC
!Read and process HVAC networks
USE MATH_FUNCTIONS, ONLY: GET_RAMP_INDEX
USE MISC_FUNCTIONS, ONLY: SEARCH_CONTROLLER
INTEGER , PARAMETER :: MAX_DUCTS = 10
INTEGER :: I1,I2,J1,J2,K1,K2,IOR
INTEGER :: IOS,IZERO,N_HVAC_READ,N,ND,ND2,NM,NN,NF,NV,I_DUCT=0,I_DUCTNODE=0,I_FAN=0 !, &
!            AIRCOIL_INDEX,I_AIRCOIL=0,FAN_INDEX, INDEX,I_FILTER=0
REAL(EB) :: AREA,DIAMETER,XYZ(3),LOSS(MAX_DUCTS,MAX_DUCTS),VOLUME_FLOW,MAX_FLOW,MAX_PRESSURE,ROUGHNESS,LENGTH,TNOW,RHOA
LOGICAL :: ROUND, SQUARE, DAMPER, REVERSE, AMBIENT
CHARACTER(30) :: TYPE_ID,DEVC_ID,CTRL_ID,DUCT_ID(MAX_DUCTS),NODE_ID(2),VENT_ID,TABLE_ID,FAN_ID,AIRCOIL_ID,RAMP_ID,ID!,SUBTYPE,FILTER_ID
CHARACTER(30), ALLOCATABLE, DIMENSION(:,:) :: NODE_DUCT_A,DUCT_NODE_A
CHARACTER(30), ALLOCATABLE, DIMENSION(:) :: NODE_FILTER_A,DUCT_FAN_A
TYPE(DUCTNODE_TYPE), POINTER :: DN=>NULL()
TYPE(DUCT_TYPE), POINTER :: DU=>NULL()
!TYPE(FILTER_TYPE), POINTER :: FI=>NULL()
!TYPE(AIRCOIL_TYPE), POINTER :: AC=>NULL()
NAMELIST /HVAC/ TYPE_ID,AREA,DIAMETER,ROUND,SQUARE,LENGTH,DAMPER,REVERSE,LOSS,DUCT_ID,NODE_ID,FAN_ID, &!FILTER_ID,AIRCOIL_ID, &
                VENT_ID,DEVC_ID,CTRL_ID,RAMP_ID,TABLE_ID,ID,VOLUME_FLOW,MAX_FLOW,MAX_PRESSURE,ROUGHNESS,XYZ,AMBIENT

TNOW=SECOND()

N_HVAC_READ = 0
 
REWIND(LU_INPUT)
COUNT_HVAC_LOOP: DO
   CALL CHECKREAD('HVAC',LU_INPUT,IOS)
   IF (IOS==1) EXIT COUNT_HVAC_LOOP
   READ(LU_INPUT,HVAC,END=15,ERR=16,IOSTAT=IOS)
   N_HVAC_READ = N_HVAC_READ + 1
   16 IF (IOS>0) THEN
         WRITE(MESSAGE,'(A,I2)') 'ERROR: Problem with HVAC line number ',N_HVAC_READ+1
         CALL SHUTDOWN(MESSAGE)
      ENDIF
   SELECT CASE(TYPE_ID)
      CASE ('DUCT')
         N_DUCTS = N_DUCTS + 1
      CASE ('NODE')
         N_DUCTNODES = N_DUCTNODES + 1
      CASE ('FILTER')
         N_FILTERS = N_FILTERS + 1
      CASE ('FAN')
         N_FANS = N_FANS + 1
      CASE ('AIRCOIL')
         N_AIRCOILS = N_AIRCOILS + 1
   END SELECT
ENDDO COUNT_HVAC_LOOP
15 CONTINUE

IF (N_HVAC_READ > 0) HVAC_SOLVE = .TRUE.

IF ((N_DUCTS > 0 .AND. N_DUCTNODES <= 0) .OR. (N_DUCTS <= 0 .AND. N_DUCTNODES > 0)) THEN
   WRITE(MESSAGE,'(A)') 'ERROR: Must have both DUCTs and DUCTNODEs in the input file'
   CALL SHUTDOWN(MESSAGE)
ENDIF

ALLOCATE(DUCT(N_DUCTS),STAT=IZERO)
CALL ChkMemErr('HVAC','DUCT',IZERO)
ALLOCATE(DUCT_NE(N_DUCTS),STAT=IZERO)
CALL ChkMemErr('HVAC','DUCT_NE',IZERO)
DUCT_NE = 0
ALLOCATE(DUCTNODE(N_DUCTNODES),STAT=IZERO)
CALL ChkMemErr('HVAC','DUCTNODE',IZERO)
ALLOCATE(DUCTNODE_NE(N_DUCTNODES),STAT=IZERO)
CALL ChkMemErr('HVAC','DUCTNODE_NE',IZERO)
DUCTNODE_NE = 0
ALLOCATE(FILTER(N_FILTERS),STAT=IZERO)
CALL ChkMemErr('HVAC','FILTER',IZERO)
ALLOCATE(FAN(N_FANS),STAT=IZERO)
CALL ChkMemErr('HVAC','FAN',IZERO)
ALLOCATE(AIRCOIL(N_AIRCOILS),STAT=IZERO)
CALL ChkMemErr('HVAC','AIRCOIL',IZERO)

!Temp arrays for input processing
ALLOCATE(DUCT_NODE_A(N_DUCTNODES,2),STAT=IZERO)
CALL ChkMemErr('HVAC','NODE_DUCTS',IZERO)
ALLOCATE(NODE_DUCT_A(N_DUCTNODES,10),STAT=IZERO)
CALL ChkMemErr('HVAC','NODE_DUCTS',IZERO)
ALLOCATE(NODE_FILTER_A(N_DUCTNODES),STAT=IZERO)
CALL ChkMemErr('HVAC','NODE_DUCTS',IZERO)
ALLOCATE(DUCT_FAN_A(N_DUCTS),STAT=IZERO)
CALL ChkMemErr('HVAC','DUCT_FAN',IZERO)

REWIND(LU_INPUT)
DO NN=1,N_HVAC_READ
   CALL SET_HVAC_DEFAULTS
   READ(LU_INPUT,HVAC)   
   SELECT CASE (TYPE_ID)
      CASE('DUCT')
         I_DUCT = I_DUCT + 1
         DU=> DUCT(I_DUCT)
         DU%ID   = ID
         IF (AREA < 0._EB) THEN
            IF (DIAMETER < 0._EB) THEN
              WRITE(MESSAGE,'(A,I2)') 'ERROR: Duct must have area or diameter specifiied, HVAC line number ',NN
              CALL SHUTDOWN(MESSAGE)
            ENDIF
            AREA = 0.5_EB*PIO2*DIAMETER**2
         ENDIF
         IF (DIAMETER < 0._EB) THEN
            IF (AREA < 0._EB) THEN
              WRITE(MESSAGE,'(A,I2)') 'ERROR: Duct must have area or diameter specified, HVAC line number ',NN
              CALL SHUTDOWN(MESSAGE)
            ENDIF
            DIAMETER = SQRT(2._EB*AREA/PIO2)
         ENDIF
         DU%AREA_INITIAL = AREA
         DU%AREA = AREA
         DU%DIAMETER = DIAMETER
         DU%LENGTH = LENGTH         
         DU%REVERSE = REVERSE
         IF (N_SPECIES > 0) THEN
            ALLOCATE(DU%YY(N_SPECIES))
            DU%YY(1:N_SPECIES) = SPECIES(1:N_SPECIES)%YY0
         ENDIF
!         IF (LOSS(1) < 0._EB .AND. VOLUME_FLOW > 1.E6_EB) THEN
!            WRITE(MESSAGE,'(A,I2)') 'ERROR: Duct needs a loss coefficient, HVAC line number ',NN
!            CALL SHUTDOWN(MESSAGE)
!         ENDIF
         DU%LOSS(1:2) = MAX(0._EB,LOSS(1:2,1))
         IF (CTRL_ID /='null' .AND. DEVC_ID /='null') THEN
            WRITE(MESSAGE,'(A,I2)') 'ERROR: Can only specify one of CTRL_ID or DEVC_ID, HVAC line number ',NN
            CALL SHUTDOWN(MESSAGE)
         ENDIF
         IF (DAMPER .AND. (FAN_ID /='null' .OR. AIRCOIL_ID /='null') .OR. &
             FAN_ID/='null' .AND. (DAMPER .OR. AIRCOIL_ID /='null') .OR. &
             AIRCOIL_ID/='null' .AND. (DAMPER .OR. FAN_ID /='null')) THEN
            WRITE(MESSAGE,'(A,I2)') 'ERROR: Duct can only have one of damper, fan or aircoil, HVAC line number ',NN
            CALL SHUTDOWN(MESSAGE)
         ENDIF
         IF (FAN_ID/='null' .AND. N_FANS<=0) THEN
            WRITE(MESSAGE,'(A,I2)') 'ERROR: Duct has fan specied but no fans have been defined, HVAC line number ',NN
            CALL SHUTDOWN(MESSAGE)
         ENDIF
         DU%DAMPER = DAMPER
         DUCT_FAN_A(I_DUCT) = FAN_ID
         IF (CTRL_ID /= 'null' .OR. DEVC_ID /= 'null') CALL SEARCH_CONTROLLER('HVAC',CTRL_ID,DEVC_ID,DU%DEVC_INDEX,DU%CTRL_INDEX,N)
         IF (DAMPER) THEN
            IF (DU%DEVC_INDEX > 0) THEN
                DU%DAMPER_OPEN = DEVICE(DU%DEVC_INDEX)%INITIAL_STATE                
            ELSEIF (DU%CTRL_INDEX > 0) THEN
                DU%DAMPER_OPEN = CONTROL(DU%CTRL_INDEX)%INITIAL_STATE                
            ELSE
                DU%DAMPER_OPEN = .TRUE.
            ENDIF
            IF (.NOT. DU%DAMPER_OPEN) DU%AREA = 0._EB
         ELSEIF (FAN_ID /='null') THEN
            IF (DU%DEVC_INDEX > 0) THEN
                DU%FAN_OPERATING = DEVICE(DU%DEVC_INDEX)%INITIAL_STATE                
            ELSEIF (DU%CTRL_INDEX > 0) THEN
                DU%FAN_OPERATING = CONTROL(DU%CTRL_INDEX)%INITIAL_STATE                
            ELSE
                DU%FAN_OPERATING = .TRUE.
            ENDIF
         ELSEIF (AIRCOIL_ID /='null') THEN
         ENDIF     
         IF (SQUARE) ROUND = .FALSE.
         DU%ROUND = ROUND
         DU%SQUARE = SQUARE
         DUCT_NODE_A(I_DUCT,:) = NODE_ID
         DU%VOLUME_FLOW_INITIAL = VOLUME_FLOW
         DU%ROUGHNESS = ROUGHNESS
      CASE('NODE')
         I_DUCTNODE = I_DUCTNODE + 1
         NODE_DUCT_A(I_DUCTNODE,:) = DUCT_ID         
         DN => DUCTNODE(I_DUCTNODE)
         DN%ID = ID
         DN%VENT_ID = VENT_ID
         DN%TABLE_ID = TABLE_ID
         DN%XYZ      = XYZ
         DN%AMBIENT  = AMBIENT
         IF (N_SPECIES > 0) THEN
            ALLOCATE(DN%YY(N_SPECIES))
            DN%YY(1:N_SPECIES) = SPECIES(1:N_SPECIES)%YY0
            ALLOCATE(DN%YY_V(N_SPECIES))            
            DN%YY_V(1:N_SPECIES) = SPECIES(1:N_SPECIES)%YY0
         ENDIF
         DO ND = 1, MAX_DUCTS
            IF (NODE_DUCT_A(I_DUCTNODE,ND) == 'null') EXIT
            DN%N_DUCTS=ND
         ENDDO
         ALLOCATE(DN%LOSS_ARRAY(MAX(2,DN%N_DUCTS),MAX(2,DN%N_DUCTS)))
         DN%LOSS_ARRAY = 0._EB
         IF (DN%N_DUCTS >=2) THEN
            DN%LOSS_ARRAY = LOSS(1:DN%N_DUCTS,1:DN%N_DUCTS)
         ELSE
            DN%LOSS_ARRAY(1,2) = LOSS(1,1)
            DN%LOSS_ARRAY(2,1) = LOSS(2,1)
         ENDIF
      CASE('FAN')
         I_FAN = I_FAN + 1
         FAN(I_FAN)%ID = ID
         FAN(I_FAN)%OFF_LOSS = LOSS(1,1)
         FAN(I_FAN)%FAN_RAMP = RAMP_ID
         FAN(I_FAN)%VOL_FLOW = VOLUME_FLOW
         FAN(I_FAN)%MAX_PRES = MAX_PRESSURE
         FAN(I_FAN)%MAX_FLOW = MAX_FLOW
         IF (RAMP_ID /= 'null') CALL GET_RAMP_INDEX(RAMP_ID,'FAN',FAN(I_FAN)%RAMP_INDEX)
         IF(( (MAX_FLOW<1.E6_EB .OR. MAX_PRESSURE<1.E6_EB) .AND. (VOLUME_FLOW<1.E6_EB .OR. RAMP_ID/='null')))THEN !.OR. &
            !(RAMP_ID/='null'.AND. (VOLUME_FLOW<1.E6_EB .OR. (MAX_FLOW<1.E6_EB .OR. MAX_PRESSURE<1.E6_EB)) .OR. &
            !(VOLUME_FLOW<1.E6_EB .AND. ((MAX_FLOW<1.E6_EB .OR. MAX_PRESSURE<1.E6_EB) .OR. RAMP_ID/='null'))) THEN
            WRITE(MESSAGE,'(A,I2)') 'ERROR: FAN can only be one of constant volume, quadratic or ramp, HVAC line number ',NN
            CALL SHUTDOWN(MESSAGE)
         ENDIF
         IF ((MAX_PRESSURE<1.E6_EB .AND. MAX_FLOW>1.E6_EB) .OR. (MAX_PRESSURE>1.E6_EB .AND. MAX_FLOW<1.E6_EB)) THEN
            WRITE(MESSAGE,'(A,I2)') 'ERROR: IF one of MAX_PRESSURE or MAX_FLOW given, both must be specified, HVAC line number '&
                                    ,NN
            CALL SHUTDOWN(MESSAGE)
         ENDIF
         IF (VOLUME_FLOW < 1.E6_EB) THEN
            FAN(I_FAN)%FAN_TYPE = 1
         ELSEIF(RAMP_ID/='null') THEN
            FAN(I_FAN)%FAN_TYPE = 3
         ELSE
            FAN(I_FAN)%FAN_TYPE = 2            
         ENDIF
   END SELECT
ENDDO

DO ND = 1, N_DUCTS
   DU => DUCT(ND)
   DO NN = 1, N_DUCTNODES 
      IF(DUCTNODE(NN)%ID == DUCT_NODE_A(ND,1)) THEN
        DU%NODE_INDEX(1) = NN
        DU%X(1) =DUCTNODE(NN)%XYZ(1)
        DU%Y(1) =DUCTNODE(NN)%XYZ(2)
        DU%Z(1) =DUCTNODE(NN)%XYZ(3)
        CYCLE
      ENDIF
      IF(DUCTNODE(NN)%ID == DUCT_NODE_A(ND,2)) THEN
        DU%NODE_INDEX(2) = NN
        DU%X(2) =DUCTNODE(NN)%XYZ(1)
        DU%Y(2) =DUCTNODE(NN)%XYZ(2)
        DU%Z(2) =DUCTNODE(NN)%XYZ(3)
        CYCLE
      ENDIF
      IF (DU%NODE_INDEX(1) > 0 .AND. DU%NODE_INDEX(2) > 0) EXIT
   ENDDO
   IF (DU%NODE_INDEX(1) <= 0) THEN
      WRITE(MESSAGE,'(A,I2)') 'ERROR: Duct node 1 not located, DUCT number ',ND
      CALL SHUTDOWN(MESSAGE)
   ENDIF
   IF (DU%NODE_INDEX(2) <= 0) THEN
      WRITE(MESSAGE,'(A,I2)') 'ERROR: Duct node 2 not located, DUCT number ',ND
      CALL SHUTDOWN(MESSAGE)
   ENDIF
   IF (DU%NODE_INDEX(1) == DU%NODE_INDEX(2)) THEN
      WRITE(MESSAGE,'(A,I2)') 'ERROR: Duct node 1 cannot be the same as duct node 2, DUCT number ',ND
      CALL SHUTDOWN(MESSAGE)
   ENDIF
   IF (DUCT_FAN_A(ND)/='null') THEN
      DO NF = 1, N_FANS
         IF(FAN(NF)%ID == DUCT_FAN_A(ND)) THEN
            DU%FAN_INDEX = NF
            EXIT
         ENDIF
      ENDDO
      IF (DU%FAN_INDEX <= 0) THEN
         WRITE(MESSAGE,'(A,I2)') 'ERROR: Fan not located for DUCT number ',ND
         CALL SHUTDOWN(MESSAGE)
      ENDIF
      IF (FAN(DU%FAN_INDEX)%FAN_TYPE==1) THEN
         DU%VOLUME_FLOW_INITIAL = FAN(DU%FAN_INDEX)%MAX_FLOW
      ENDIF
   ENDIF   
ENDDO

RHOA = P_INF/(TMPA*RSUM0)

DO NN = 1, N_DUCTNODES
   DN => DUCTNODE(NN)
   IF (DN%VENT_ID /= 'null') THEN
      MESH_LOOP: DO NM = 1, NMESHES
         DO NV = 1, MESHES(NM)%N_VENT
            IF(MESHES(NM)%VENTS(NV)%ID == DN%VENT_ID) THEN
               IF (MESHES(NM)%VENTS(NV)%IBC /= HVAC_SURF_INDEX) THEN
                  WRITE(MESSAGE,'(A,A)') 'ERROR: DUCTNODE attached to VENT without SURF_ID HVAC for VENT ID ',&
                                          MESHES(NM)%VENTS(NV)%ID
               ENDIF
               DN%MESH_INDEX = NM
               DN%VENT_INDEX = NV
               IOR = MESHES(NM)%VENTS(NV)%IOR
               I1 = MESHES(NM)%VENTS(NV)%I1
               I2 = MESHES(NM)%VENTS(NV)%I2
               J1 = MESHES(NM)%VENTS(NV)%J1
               J2 = MESHES(NM)%VENTS(NV)%J2
               K1 = MESHES(NM)%VENTS(NV)%K1
               K2 = MESHES(NM)%VENTS(NV)%K2
               SELECT CASE (ABS(IOR))
                  CASE (1)
                     J1 = J1 + 1
                     K1 = K1 + 1
                  CASE (2)
                     I1 = I1 + 1
                     K1 = K1 + 1
                  CASE (3)
                     I1 = I1 + 1
                     J1 = J1 + 1
               END SELECT               
               MESHES(NM)%VENTS(NV)%NODE_INDEX=NN
               DN%XYZ(1) = 0.5_EB*(MESHES(NM)%X(I1)+MESHES(NM)%X(I2))
               DN%XYZ(2) = 0.5_EB*(MESHES(NM)%Y(J1)+MESHES(NM)%Y(J2))
               DN%XYZ(3) = 0.5_EB*(MESHES(NM)%Z(K1)+MESHES(NM)%Z(K2))
               EXIT MESH_LOOP            
            ENDIF
         ENDDO
      ENDDO MESH_LOOP
   ENDIF
   IF (DN%VENT_INDEX > 0 .AND. DN%AMBIENT) THEN
      WRITE(MESSAGE,'(A,I2)') 'ERROR: DUCTNODE cannot be AMBIENT and have an assigned VENT_ID, HVAC line number ',N
      CALL SHUTDOWN(MESSAGE)
   ENDIF
   IF (DN%N_DUCTS == 1 .AND. DN%VENT_INDEX < 0 .AND. .NOT. DN%AMBIENT) THEN
      WRITE(MESSAGE,'(A,I2)') 'ERROR: Internal DUCTNODE must have at least two attached ducts, DUCTNODE number ',NN
      CALL SHUTDOWN(MESSAGE)
   ENDIF  
   IF (DN%N_DUCTS> 1 .AND. (DN%AMBIENT .OR. DN%VENT_INDEX > 0) ) THEN
      WRITE(MESSAGE,'(A,I2)') 'ERROR: External DUCTNODE can only have one attached duct, DUCTNODE number ',NN
      CALL SHUTDOWN(MESSAGE)
   ENDIF  
   ALLOCATE(DN%DUCT_INDEX(DN%N_DUCTS))
   ALLOCATE(DN%DIR(DN%N_DUCTS))
   DN%DUCT_INDEX = -1
   ALLOCATE(DN%BRANCH_LOSS(DN%N_DUCTS))   
   DN%BRANCH_LOSS = 0._EB
   DO ND = 1,DN%N_DUCTS
      DO ND2 = 1, N_DUCTS
         IF (NODE_DUCT_A(NN,ND) == DUCT(ND2)%ID) THEN
            DN%DUCT_INDEX(ND) = ND2
            IF (DUCT(ND2)%NODE_INDEX(1)==NN) THEN
               DN%DIR(ND) = 1._EB
            ELSE
               DN%DIR(ND) =-1._EB
            ENDIF
            CYCLE
         ENDIF
      ENDDO
      IF (DN%DUCT_INDEX(ND)==-1) THEN
         WRITE(MESSAGE,'(A,I2,A,I2)') 'ERROR: DUCT ',ND,' not found for DUCTNODE number ',NN
         CALL SHUTDOWN(MESSAGE)
      ENDIF  
   ENDDO
   
   !Initialize duct node properties
   DN%TMP    = TMPA
   DN%TMP_V  = TMPA
   DN%RSUM   = RSUM0
   DN%RSUM_V = RSUM0
   DN%CP_V   = Y2CP_C(NINT(TMPA))
   DN%CP     = Y2CP_C(NINT(TMPA))
   IF (STRATIFICATION) THEN
      DN%P     =  P_INF*EXP(GVEC(3)*(DN%XYZ(3)-GROUND_LEVEL)/(RSUM0*TMPA))
      DN%RHO   =  DN%P/(TMPA*RSUM0)
      DN%RHO_V  = DN%RHO
   ELSE
      DN%P      = P_INF
      DN%RHO    = RHOA
      DN%RHO_V  = RHOA
   ENDIF
ENDDO

!Temp arrays for input processing
IF (ALLOCATED(DUCT_NODE_A)) DEALLOCATE(DUCT_NODE_A)
IF (ALLOCATED(NODE_DUCT_A)) DEALLOCATE(NODE_DUCT_A)
IF (ALLOCATED(NODE_FILTER_A)) DEALLOCATE(NODE_FILTER_A)
IF (ALLOCATED(DUCT_FAN_A)) DEALLOCATE(DUCT_FAN_A)

CALL DETERMINE_FIXED_ELEMENTS

TUSED(16,1)=TUSED(16,1)+SECOND()-TNOW

RETURN

CONTAINS 

SUBROUTINE SET_HVAC_DEFAULTS

TYPE_ID = 'null'
TABLE_ID = 'null'
CTRL_ID = 'null'
DEVC_ID = 'null'
DUCT_ID = 'null'
NODE_ID = 'null'
VENT_ID = 'null'
RAMP_ID = 'null'
FAN_ID = 'null'
AIRCOIL_ID = 'null'
AREA = -1._EB
DIAMETER = -1._EB
AMBIENT = .FALSE.
ROUND = .TRUE.
SQUARE = .FALSE.
LENGTH = 1._EB
DAMPER = .FALSE.
REVERSE = .FALSE.
LOSS = 0._EB
VOLUME_FLOW  = 1.E7_EB
MAX_FLOW     = 1.E7_EB
MAX_PRESSURE = 1.E7_EB
ROUGHNESS    = 0._EB
XYZ          = 0.
RETURN

END SUBROUTINE SET_HVAC_DEFAULTS

END SUBROUTINE READ_HVAC

SUBROUTINE HVAC_CALC(T)
!Solve for flows in the HVAC networks
INTEGER :: NNE
REAL(EB) :: T
LOGICAL :: CHANGE=.TRUE.,DPSTARCALL
TYPE(NETWORK_TYPE), POINTER:: NE=>NULL()
DPSTARCALL = .FALSE.
CALL FIND_NETWORKS(CHANGE)
!WRITE(*,*) 'ENTER HVAC_CALC'
DO NNE = 1, N_NETWORKS
   NE =>NETWORK(NNE)
   CALL SET_GUESS(NNE,ITER)  
   IF (NE%N_MATRIX > 0) THEN 
      ALLOCATE(LHS(NE%N_MATRIX,NE%N_MATRIX))
      ALLOCATE(RHS(NE%N_MATRIX))
      IF (.NOT. DPSTARCALL .AND. N_ZONE > 0) THEN
         CALL DPSTARCALC
         DPSTARCALL = .TRUE.
      ENDIF
      ITER = 0
      DO WHILE (ITER < ITER_MAX)
!         WRITE(*,*) '**** ITER **** ',ITER
         LHS = 0._EB
         RHS = 0._EB      
         CALL SET_GUESS(NNE,ITER)      
         CALL SET_DONOR(NNE,T)
         CALL UPDATE_LOSS(NNE)
         CALL RHSNODE(NNE)
         CALL RHSDUCT(NNE)   
         CALL LHSNODE(NNE)
         CALL LHSDUCT(NNE)
         CALL MATRIX_SOLVE(NNE)
         CALL HVAC_UPDATE(NNE)
         CALL CONVERGENCE_CHECK(ITER,NNE)
         ITER = ITER + 1
      ENDDO
      DEALLOCATE(LHS)
      DEALLOCATE(RHS)
   ELSE
      CALL SET_DONOR(NNE,T)
      CALL HVAC_UPDATE(NNE)
   ENDIF
ENDDO
IF (CORRECTOR) DUCT%VEL(OLD) = DUCT%VEL(NEW)

END SUBROUTINE HVAC_CALC

SUBROUTINE MATRIX_SOLVE(NNE)
USE MATH_FUNCTIONS,ONLY : GAUSSJ
INTEGER :: NNE,IERR,ND,NN
TYPE(NETWORK_TYPE), POINTER :: NE=>NULL()
TYPE(DUCT_TYPE), POINTER :: DU=>NULL()
TYPE(DUCTNODE_TYPE), POINTER :: DN=>NULL()

NE =>NETWORK(NNE)
!WRITE(*,'(2X,6(E11.5,1X),A)') (LHS(ND,1),ND=1,NE%N_MATRIX),' LHS'
!WRITE(*,'(2X,6(E11.5,1X),A)') (LHS(ND,2),ND=1,NE%N_MATRIX),' LHS'
!WRITE(*,'(2X,6(E11.5,1X),A)') (LHS(ND,3),ND=1,NE%N_MATRIX),' LHS'
!WRITE(*,'(2X,6(E11.5,1X),A)') (LHS(ND,4),ND=1,NE%N_MATRIX),' LHS'
!WRITE(*,'(2X,6(E11.5,1X),A)') (LHS(ND,5),ND=1,NE%N_MATRIX),' LHS'
!WRITE(*,'(2X,6(E11.5,1X),A)') (LHS(ND,6),ND=1,NE%N_MATRIX),' LHS'
!WRITE(*,'(2X,6(E11.5,1X),A)') RHS,' RHS_IN'
!WRITE(*,*) RHS,' RHS_IN'
CALL GAUSSJ(LHS,NE%N_MATRIX,NE%N_MATRIX,RHS,1,1,IERR)
DO ND = 1,NE%N_DUCTS
   DU=>DUCT(NE%DUCT_INDEX(ND))
   IF (DU%FIXED .OR. DU%AREA < ZERO_P) CYCLE
   DU%VEL(NEW) = RHS(NE%MATRIX_INDEX(ND))
ENDDO
DO NN = 1,NE%N_DUCTNODES
   DN=>DUCTNODE(NE%NODE_INDEX(NN))
   IF (DN%FIXED) CYCLE
   DN%P = RHS(NE%MATRIX_INDEX(NE%N_DUCTS+NN))
ENDDO
!WRITE(*,'(2X,6(E11.5,1X),A)') RHS,' RHS_OUT'
!WRITE(*,*) RHS,' RHS_OUT'

END SUBROUTINE MATRIX_SOLVE

SUBROUTINE HVAC_UPDATE(NNE)
!Iterate duct network to update all ducts and nodes
USE COMP_FUNCTIONS, ONLY: SECOND
USE PHYSICAL_FUNCTIONS, ONLY : GET_AVERAGE_SPECIFIC_HEAT,GET_SPECIFIC_GAS_CONSTANT
REAL(EB) :: TNOW,MTOT,ETOT,VTOT,YYTOT(1:N_SPECIES),TGUESS,VFLOW,YY_N(1:N_SPECIES)
INTEGER, INTENT(IN) :: NNE
INTEGER :: NN,ND!,COUNTER1,COUNTER2
LOGICAL :: CYCLE_FLAG
TYPE (DUCTNODE_TYPE), POINTER :: DN=>NULL()
TYPE (DUCT_TYPE), POINTER :: DU=>NULL()
TYPE (NETWORK_TYPE), POINTER :: NE=>NULL()

TNOW=SECOND()
NE => NETWORK(NNE)
DUCT%UPDATED = .FALSE.
DO NN = 1,NE%N_DUCTNODES
   DN=>DUCTNODE(NE%NODE_INDEX(NN))
   DN%UPDATED = .FALSE.
   IF (DN%VENT_INDEX > 0) THEN   
      DU => DUCT(DN%DUCT_INDEX(1))
      IF (DU%VEL(NEW)*DN%DIR(1) > 0._EB) DN%UPDATED = .TRUE.
   ENDIF
ENDDO

DO ND = 1,NE%N_DUCTS
   DUCT(NE%DUCT_INDEX(ND))%UPDATED = .FALSE.
ENDDO

ITER_LOOP: DO
   CYCLE_FLAG = .FALSE.
   DUCT_LOOP:DO ND = 1,NE%N_DUCTS
      DU=>DUCT(NE%DUCT_INDEX(ND))
      IF (DU%UPDATED) CYCLE DUCT_LOOP
      CYCLE_FLAG = .TRUE.
      IF (DU%VEL(NEW) > ZERO_P) THEN
         DN => DUCTNODE(DU%NODE_INDEX(1))
      ELSEIF (DU%VEL(NEW) < -ZERO_P) THEN
         DN => DUCTNODE(DU%NODE_INDEX(2))
      ELSE
         DU%UPDATED = .TRUE.
         DU%VEL(NEW) = 0._EB
         CYCLE DUCT_LOOP
      ENDIF
      IF (DN%UPDATED) THEN
         DU%RHO_D  = DN%RHO
         DU%TMP_D  = DN%TMP
         DU%CP_D   = DN%CP         
         DU%UPDATED = .TRUE.
         IF (N_SPECIES > 0) DU%YY(:) = DN%YY(:)
      ENDIF
   ENDDO DUCT_LOOP
   NODE_LOOP:DO NN = 1,NE%N_DUCTNODES
      DN=>DUCTNODE(NE%NODE_INDEX(NN))
      IF(DN%UPDATED) CYCLE NODE_LOOP
      CYCLE_FLAG = .TRUE.
      VTOT = 0._EB
      MTOT = 0._EB
      ETOT = 0._EB
      IF(N_SPECIES > 0) YYTOT = 0._EB
      TGUESS = 0._EB
      DO ND = 1,DN%N_DUCTS
         DU => DUCT(DN%DUCT_INDEX(ND))         
         IF (DU%AREA<ZERO_P) CYCLE
         IF (DU%VEL(NEW)*DN%DIR(ND) > 0._EB) CYCLE
         IF (DU%VEL(NEW)*DN%DIR(ND) < 0._EB .AND. .NOT. DU%UPDATED) CYCLE NODE_LOOP
         VFLOW = ABS(DU%VEL(NEW)*DU%AREA)
         VTOT = VTOT + VFLOW
         MTOT = MTOT + VFLOW * DU%RHO_D
         ETOT = ETOT + VFLOW * DU%RHO_D * DU%TMP_D * DU%CP_D
         IF(N_SPECIES > 0) YYTOT = YYTOT + VFLOW * DU%RHO_D * DU%YY
         TGUESS = TGUESS + VFLOW * DU%RHO_D * DU%TMP_D
      ENDDO
      DN%UPDATED = .TRUE.
      IF (ABS(MTOT)<=ZERO_P) CYCLE NODE_LOOP
      IF(N_SPECIES > 0) THEN
         DN%YY(:)  = YYTOT/MTOT
         YY_N = DN%YY
         CALL GET_SPECIFIC_GAS_CONSTANT(YY_N,DN%RSUM)
      ENDIF
      TGUESS = TGUESS / MTOT      
      ETOT = ETOT/ MTOT
      DN%TMP = TGUESS
      CP_LOOP: DO
         IF (N_SPECIES>0) THEN
            CALL GET_AVERAGE_SPECIFIC_HEAT(YY_N,DN%CP,MIN(5000,NINT(TGUESS)))
         ELSE
            DN%CP = Y2CP_C(MIN(5000,NINT(TGUESS)))
         ENDIF
         DN%TMP = ETOT /DN%CP
         IF (ABS(DN%TMP - TGUESS) < ZERO_P) EXIT CP_LOOP
         IF ((DN%TMP - TGUESS)/DN%TMP < 0.0005_EB) EXIT CP_LOOP
         TGUESS = DN%TMP
      ENDDO CP_LOOP
      DN%RHO = DN%P/(DN%RSUM*DN%TMP)
   ENDDO NODE_LOOP
   IF (.NOT. CYCLE_FLAG) EXIT ITER_LOOP
ENDDO ITER_LOOP
!WRITE(*,*) DUCT%VEL(NEW),' V'
!WRITE(*,*) DUCTNODE%TMP,' T'
!WRITE(*,*) DUCTNODE%P,' P'
!WRITE(*,*) (DUCTNODE%RHO),' R'
!WRITE(*,*) (DUCTNODE%RSUM),' RSUM'
!IF (N_SPECIES > 0) WRITE(*,*) (DUCTNODE(NN)%ID,DUCTNODE(NN)%YY(2),NN=1,N_DUCTNODES)

!WRITE(*,*) 'EXIT CALC'

TUSED(16,:)=TUSED(16,:)+SECOND()-TNOW

END SUBROUTINE HVAC_UPDATE


SUBROUTINE RHSNODE(NETWORK_INDEX)
USE GLOBAL_CONSTANTS
INTEGER, INTENT(IN)::NETWORK_INDEX
INTEGER :: NN,ND, ARRAYLOC
TYPE(NETWORK_TYPE), POINTER::NE=>NULL()
TYPE(DUCT_TYPE), POINTER::DU=>NULL()
TYPE(DUCTNODE_TYPE), POINTER::DN=>NULL()

NE => NETWORK(NETWORK_INDEX)
DO NN = 1, NE%N_DUCTNODES
   DN => DUCTNODE(NE%NODE_INDEX(NN))
   IF (DN%FIXED .OR. DN%VENT_INDEX > 0 .OR. DN%AMBIENT) CYCLE
   ARRAYLOC = NE%MATRIX_INDEX(NE%N_DUCTS+DUCTNODE_NE(NE%NODE_INDEX(NN)))
   DO ND = 1,DN%N_DUCTS
      DU => DUCT(DN%DUCT_INDEX(ND))
      IF (DU%FIXED) RHS(ARRAYLOC) = RHS(ARRAYLOC) + DN%DIR(ND)*DU%RHO_D*DU%VOLUME_FLOW
   END DO
ENDDO

END SUBROUTINE RHSNODE

SUBROUTINE LHSNODE(NETWORK_INDEX)
USE GLOBAL_CONSTANTS
INTEGER, INTENT(IN)::NETWORK_INDEX
INTEGER :: NN,ND, ARRAYLOC1,ARRAYLOC2
TYPE(NETWORK_TYPE), POINTER::NE=>NULL()
TYPE(DUCT_TYPE), POINTER::DU=>NULL()
TYPE(DUCTNODE_TYPE), POINTER::DN=>NULL()

NE => NETWORK(NETWORK_INDEX)
DO NN = 1, NE%N_DUCTNODES
   DN => DUCTNODE(NE%NODE_INDEX(NN))
   IF (DN%FIXED .OR. DN%VENT_INDEX > 0) CYCLE
   ARRAYLOC1 = NE%MATRIX_INDEX(NE%N_DUCTS+DUCTNODE_NE(NE%NODE_INDEX(NN)))
   DO ND = 1,DN%N_DUCTS
      DU => DUCT(DN%DUCT_INDEX(ND))
      IF (DU%FIXED .OR. DU%AREA <ZERO_P) CYCLE
      ARRAYLOC2 = NE%MATRIX_INDEX(DUCT_NE(DN%DUCT_INDEX(ND)))
      LHS(ARRAYLOC1,ARRAYLOC2) = -DN%DIR(ND)*DU%RHO_D*DU%AREA      
   END DO
ENDDO

END SUBROUTINE LHSNODE

SUBROUTINE DPSTARCALC
USE GLOBAL_CONSTANTS
INTEGER :: NN,IPZ,NM
TYPE(DUCT_TYPE), POINTER::DU=>NULL()
TYPE(DUCTNODE_TYPE), POINTER::DN=>NULL()
TYPE(P_ZONE_TYPE), POINTER::PZ=>NULL()

DO IPZ = 1,N_ZONE
   PZ => P_ZONE(IPZ)
   IF (PZ%N_DUCTNODES==0) CYCLE
   P_ZONE(IPZ)%DPSTAR = P_ZONE(IPZ)%DPSTAR * DT
   DO NN = 1,PZ%N_DUCTNODES
      DN=>DUCTNODE(PZ%NODE_INDEX(NN))
      NM = DN%MESH_INDEX
      DU=>DUCT(DN%DUCT_INDEX(1))
      P_ZONE(IPZ)%DPSTAR = P_ZONE(IPZ)%DPSTAR  + DN%DIR(1) * DU%AREA * DU%VEL(OLD) * DT/PSUM(IPZ,NM)
      IF (DU%FIXED) P_ZONE(IPZ)%DPSTAR = P_ZONE(IPZ)%DPSTAR  - DN%DIR(1) * DU%AREA * DU%VEL(NEW) * DT/PSUM(IPZ,NM)
   ENDDO
ENDDO

END SUBROUTINE DPSTARCALC

SUBROUTINE RHSDUCT(NETWORK_INDEX)
USE GLOBAL_CONSTANTS
INTEGER, INTENT(IN)::NETWORK_INDEX
INTEGER :: ND, ARRAYLOC,IPZ,NM
REAL(EB) :: HEAD,XYZ(3)
TYPE(NETWORK_TYPE), POINTER::NE=>NULL()
TYPE(DUCT_TYPE), POINTER::DU=>NULL()
TYPE(DUCTNODE_TYPE), POINTER::DN=>NULL()
TYPE(P_ZONE_TYPE), POINTER::PZ=>NULL()

NE => NETWORK(NETWORK_INDEX)
DO ND = 1, NE%N_DUCTS
   DU => DUCT(NE%DUCT_INDEX(ND))
   IF (DU%FIXED .OR. DU%AREA < ZERO_P) CYCLE
   HEAD = 0._EB
   ARRAYLOC = NE%MATRIX_INDEX(DUCT_NE(NE%DUCT_INDEX(ND)))
   IF (DU%FAN_INDEX>0) CALL UPDATE_FAN(ND)
   DN=>DUCTNODE(DU%NODE_INDEX(1))
   IF (DN%AMBIENT) THEN
      HEAD = HEAD + DN%P
   ELSEIF (DN%VENT_INDEX > 0) THEN
      HEAD = HEAD + DN%P   
!      WRITE(*,*) 'HEAD  ',HEAD  ,DN%P             
      IF (N_ZONE > 0) THEN
         IPZ = DN%ZONE_INDEX
         NM = DN%MESH_INDEX
         IF (IPZ > 0) THEN
            PZ => P_ZONE(IPZ)
            HEAD = HEAD + PZ%DPSTAR
!            WRITE(*,*) 'HEAD1 ',HEAD , PZ%DPSTAR
         ENDIF
      ENDIF
   ENDIF
   XYZ = DN%XYZ
   DN=>DUCTNODE(DU%NODE_INDEX(2))
   IF (DN%AMBIENT) THEN
      HEAD = HEAD - DN%P
   ELSEIF (DN%VENT_INDEX > 0) THEN
      HEAD = HEAD - DN%P   
!      WRITE(*,*) 'HEADA  ',HEAD  ,DN%P             
      IF (N_ZONE > 0) THEN
         IPZ = DN%ZONE_INDEX
         NM = DN%MESH_INDEX
         IF (IPZ > 0) THEN
            PZ => P_ZONE(IPZ)
            HEAD = HEAD - PZ%DPSTAR
!            WRITE(*,*) 'HEAD1A ',HEAD , PZ%DPSTAR
         ENDIF
      ENDIF
   ENDIF
   XYZ = DN%XYZ - XYZ
   IF (STRATIFICATION) THEN
      HEAD = HEAD + DN%P*(1._EB-EXP(-(GVEC(1)*XYZ(1)+GVEC(2)*XYZ(2)+GVEC(3)*XYZ(3))/(DU%RSUM_D*DU%TMP_D)))
!      WRITE(*,*) 'HEAD5-',XYZ,DU%RSUM_D
!      WRITE(*,*) 'HEAD5 ',HEAD,DN%P*(1._EB-EXP(-(GVEC(1)*XYZ(1)+GVEC(2)*XYZ(2)+GVEC(3)*XYZ(3))/(DU%RSUM_D*DU%TMP_D))),'A'
   ELSE
      HEAD = HEAD + (GVEC(1)*XYZ(1)+GVEC(2)*XYZ(2)+GVEC(3)*XYZ(3))*DU%RHO_D
!      WRITE(*,*) 'HEAD5 ',HEAD,(GVEC(1)*XYZ(1)+GVEC(2)*XYZ(2)+GVEC(3)*XYZ(3))*DU%RHO_D
   ENDIF
   RHS(ARRAYLOC) = DU%VEL(OLD)+DT/DU%LENGTH*((HEAD+DU%DP_FAN)/DU%RHO_D-0.5_EB*DU%TOTAL_LOSS*ABS(DU%VEL(NEW))*DU%VEL(GUESS))
   !WRITE(*,*) 'RHSD ',ND,DU%VEL(OLD),HEAD,DU%DP_FAN,0.5_EB*DU%TOTAL_LOSS*ABS(DU%VEL(NEW))*DU%VEL(GUESS)
!   DN=>DUCTNODE(DU%NODE_INDEX(2))
!   IF (DN%VENT_INDEX > 0) THEN
!      IF (DN%ZONE_INDEX >0) THEN
!         PZ => P_ZONE(DN%ZONE_INDEX)
!         DO NN=1,PZ%N_DUCTNODES
!            DN2=>DUCTNODE(PZ%NODE_INDEX(NN))
!            DU2=>DUCT(DN%DUCT_INDEX(1))
!            ARRAYLOC2 = NE%MATRIX_INDEX(DUCT_NE(DN2%DUCT_INDEX(1)))
!            LHS(ARRAYLOC1,ARRAYLOC2) = LHS(ARRAYLOC1,ARRAYLOC2) - PZ%TMP/PZ%VOLUME * DN2%DIR(1) * DU2%RHO_D * DU2%RSUM_D * DU2%AREA
!         ENDDO
!      ENDIF
!   ENDIF
   
ENDDO

END SUBROUTINE RHSDUCT


SUBROUTINE LHSDUCT(NETWORK_INDEX)
USE GLOBAL_CONSTANTS
INTEGER, INTENT(IN)::NETWORK_INDEX
INTEGER :: NN,ND, ARRAYLOC1,ARRAYLOC2
TYPE(NETWORK_TYPE), POINTER::NE=>NULL()
TYPE(DUCT_TYPE), POINTER::DU=>NULL(),DU2=>NULL()
TYPE(DUCTNODE_TYPE), POINTER::DN=>NULL(),DN2=>NULL()
TYPE(P_ZONE_TYPE), POINTER::PZ=>NULL()
NE => NETWORK(NETWORK_INDEX)

DO ND = 1, NE%N_DUCTS
   DU => DUCT(NE%DUCT_INDEX(ND))
   IF (DU%FIXED .OR. DU%AREA < ZERO_P) CYCLE
   ARRAYLOC1 = NE%MATRIX_INDEX(DUCT_NE(NE%DUCT_INDEX(ND)))
   LHS(ARRAYLOC1,ARRAYLOC1) = 1._EB+0.5_EB*DU%TOTAL_LOSS*ABS(DU%VEL(NEW)-DU%VEL(GUESS))*DT/DU%LENGTH
   !LHS(ARRAYLOC1,ARRAYLOC1) = 1._EB+0.5_EB*DU%TOTAL_LOSS*ABS(DU%VEL(NEW))*DT/DU%LENGTH
   DN=>DUCTNODE(DU%NODE_INDEX(1))
   IF (DN%VENT_INDEX < 0) THEN
      IF (.NOT. DN%AMBIENT) THEN
         ARRAYLOC2 = NE%MATRIX_INDEX(NE%N_DUCTS+DUCTNODE_NE(DU%NODE_INDEX(1)))
         LHS(ARRAYLOC1,ARRAYLOC2) = -DT/(DU%RHO_D*DU%LENGTH)
      ENDIF
   ELSE
      IF (DN%ZONE_INDEX >0) THEN
         PZ => P_ZONE(DN%ZONE_INDEX)
         DO NN=1,PZ%N_DUCTNODES
            DN2=>DUCTNODE(PZ%NODE_INDEX(NN))
            DU2=>DUCT(DN2%DUCT_INDEX(1))
            IF (DU2%AREA < ZERO_P .OR. DU2%FIXED) CYCLE
            ARRAYLOC2 = NE%MATRIX_INDEX(DUCT_NE(DN2%DUCT_INDEX(1)))
            LHS(ARRAYLOC1,ARRAYLOC2) = LHS(ARRAYLOC1,ARRAYLOC2) + DN2%DIR(1)*DU2%AREA*DT**2 / &
                                       (PSUM(DN%ZONE_INDEX,DN%MESH_INDEX)*DU%RHO_D*DU%LENGTH)
         ENDDO
      ENDIF
   ENDIF
   DN=>DUCTNODE(DU%NODE_INDEX(2))
   IF (DN%VENT_INDEX < 0) THEN
      IF (.NOT. DN%AMBIENT) THEN
         ARRAYLOC2 = NE%MATRIX_INDEX(NE%N_DUCTS+DUCTNODE_NE(DU%NODE_INDEX(2)))         
         LHS(ARRAYLOC1,ARRAYLOC2) = DT/(DU%RHO_D*DU%LENGTH)
      ENDIF
   ELSE
      IF (DN%ZONE_INDEX >0) THEN
         PZ => P_ZONE(DN%ZONE_INDEX)
         DO NN=1,PZ%N_DUCTNODES
            DN2=>DUCTNODE(PZ%NODE_INDEX(NN))
            DU2=>DUCT(DN2%DUCT_INDEX(1))
            IF (DU2%AREA < ZERO_P .OR. DU2%FIXED) CYCLE  
            ARRAYLOC2 = NE%MATRIX_INDEX(DUCT_NE(DN2%DUCT_INDEX(1)))
            LHS(ARRAYLOC1,ARRAYLOC2) = LHS(ARRAYLOC1,ARRAYLOC2) - DN2%DIR(1)*DU2%AREA*DT**2 / &
                                       (PSUM(DN%ZONE_INDEX,DN%MESH_INDEX)*DU%RHO_D*DU%LENGTH)
         ENDDO
      ENDIF
   ENDIF
ENDDO

END SUBROUTINE LHSDUCT

SUBROUTINE UPDATE_FAN(DUCT_INDEX)
USE MATH_FUNCTIONS, ONLY : EVALUATE_RAMP
INTEGER, INTENT(IN) :: DUCT_INDEX
!TYPE (DUCTNODE_TYPE), POINTER :: DN=>NULL()
TYPE(DUCT_TYPE), POINTER::DU=>NULL()
TYPE(FAN_TYPE), POINTER::FA=>NULL()
REAL(EB) :: VOLUME_FLUX,DP_OLD

DU=> DUCT(DUCT_INDEX)
FA=> FAN(DU%FAN_INDEX)

IF (DU%DEVC_INDEX > 0) THEN
   DU%FAN_OPERATING = DEVICE(DU%DEVC_INDEX)%CURRENT_STATE
ELSEIF (DU%CTRL_INDEX > 0) THEN
   DU%FAN_OPERATING = CONTROL(DU%CTRL_INDEX)%CURRENT_STATE
ENDIF

IF (.NOT. DU%FAN_OPERATING) THEN
   DU%DP_FAN = 0._EB
   RETURN
ENDIF
DP_OLD = DU%DP_FAN
SELECT CASE (FA%FAN_TYPE)
   CASE(1) !Constant flow
      RETURN
   CASE(2) !Quadratic
      VOLUME_FLUX = DU%VEL(GUESS)*DU%AREA
      IF (DU%REVERSE) VOLUME_FLUX = -VOLUME_FLUX      
      IF (FA%MAX_FLOW < 0._EB) THEN
         IF (VOLUME_FLUX > 0._EB) THEN
            DU%DP_FAN = -FA%MAX_PRES
         ELSE
            DU%DP_FAN = -MAX(0._EB,FA%MAX_PRES * (1._EB - (VOLUME_FLUX/FA%MAX_FLOW)**2))
         ENDIF
      ELSE
         IF (VOLUME_FLUX < 0._EB) THEN
            DU%DP_FAN = FA%MAX_PRES
         ELSE
            DU%DP_FAN = MAX(0._EB,FA%MAX_PRES * (1._EB - (VOLUME_FLUX/FA%MAX_FLOW)**2))        
         ENDIF
      ENDIF
   CASE(3) !Fan curve
      VOLUME_FLUX = DU%VEL(GUESS)*DU%AREA
      IF (DU%REVERSE) VOLUME_FLUX = -VOLUME_FLUX
      DU%DP_FAN = EVALUATE_RAMP(VOLUME_FLUX,0._EB,FA%RAMP_INDEX)
END SELECT
IF (DU%REVERSE) DU%DP_FAN = -DU%DP_FAN
DU%DP_FAN = 0.3_EB*DU%DP_FAN+0.7_EB*DP_OLD

END SUBROUTINE UPDATE_FAN


SUBROUTINE HVAC_BC_IN(NM)
!Average gas properties at VENTs connected to HVAC system
USE PHYSICAL_FUNCTIONS, ONLY : GET_SPECIFIC_GAS_CONSTANT,GET_AVERAGE_SPECIFIC_HEAT,GET_SPECIFIC_ENTHALPY
INTEGER, INTENT(IN) :: NM
INTEGER :: NN,II,JJ,KK,I1,I2,J1,J2,K1,K2,IC,IW,IOR,NODE_ZONE
REAL(EB) :: YY_SUM(N_SPECIES),YY_N(1:N_SPECIES),RHO_SUM,P_SUM,AREA_SUM,AREA,TMP_SUM,H_SUM,H_G,TMP_NEW,CP,TNOW,VEL2,P_AVE
REAL(EB), POINTER, DIMENSION(:,:,:) :: RHOP,UP,VP,WP
REAL(EB), POINTER, DIMENSION(:,:) :: PBARP
TYPE (MESH_TYPE),POINTER :: M=>NULL()
TYPE (VENTS_TYPE),POINTER :: VT=>NULL()
TYPE (DUCTNODE_TYPE), POINTER :: DN=>NULL()

IF (EVACUATION_ONLY(NM)) RETURN

TNOW=SECOND()

M => MESHES(NM)
CALL POINT_TO_MESH(NM)

IF (PREDICTOR) THEN
   PBARP => PBAR
   RHOP  => RHO
   UP    => U
   VP    => V
   WP    => W      
ELSE
   PBARP => PBAR_S
   RHOP  => RHOS
   UP    => US
   VP    => VS
   WP    => WS      
ENDIF

DO NN = 1, N_DUCTNODES
   IF (N_SPECIES>0) YY_SUM=0._EB
   RHO_SUM=0._EB
   P_SUM=0._EB
   AREA_SUM=0._EB
   TMP_SUM=0._EB
   H_SUM=0._EB
   NODE_ZONE=-1
   DN=>DUCTNODE(NN)
   IF (DN%MESH_INDEX /= NM .OR. DN%VENT_INDEX < 0) CYCLE
   IF (NM==1 .AND. DN%AMBIENT) THEN
      IF (STRATIFICATION) THEN
         DN%TMP = TMPA + LAPSE_RATE*DN%XYZ(3)
         IF (ABS(LAPSE_RATE) > ZERO_P) THEN
            DN%P = P_INF*(DN%TMP/TMPA)**(GVEC(3)/RSUM0/LAPSE_RATE)
         ELSE
            DN%P = P_INF*EXP(GVEC(3)*DN%XYZ(3)/(RSUM0*TMPA))
         ENDIF
      ELSE
         DN%TMP = TMPA
         DN%P   = P_INF
      ENDIF
      IF (N_SPECIES > 0) THEN
         DN%YY(1:N_SPECIES) = SPECIES(1:N_SPECIES)%YY0
         YY_N = DN%YY
         CALL GET_AVERAGE_SPECIFIC_HEAT(YY_N,DN%CP,MIN(5000,NINT(DN%TMP)))
         CALL GET_SPECIFIC_GAS_CONSTANT(YY_N,DN%RSUM_V)
      ELSE
         DN%CP_V = Y2CP_C(NINT(DN%TMP))
         DN%RSUM_V = RSUM0
      ENDIF
      DN%RHO = DN%P/(DN%TMP*DN%RSUM_V)
   ENDIF
   VT=>VENTS(DN%VENT_INDEX)
   IOR = VT%IOR
   I1 = VT%I1
   I2 = VT%I2
   J1 = VT%J1
   J2 = VT%J2
   K1 = VT%K1
   K2 = VT%K2
   SELECT CASE (ABS(IOR))
      CASE (1)
         J1 = J1 + 1
         K1 = K1 + 1
         IF (IOR > 0) THEN
            I1 = I1 + 1
            I2 = I1
         ENDIF
      CASE (2)
         I1 = I1 + 1
         K1 = K1 + 1
         IF (IOR > 0) THEN
            J1 = J1 + 1
            J2 = J1
         ENDIF
      CASE (3)
         I1 = I1 + 1
         J1 = J1 + 1
         IF (IOR > 0) THEN
            K1 = K1 + 1
            K2 = K1
         ENDIF
   END SELECT
   DO KK = K1,K2
      DO JJ = J1,J2
         DO II = I1, I2
            IC = CELL_INDEX(II,JJ,KK)
            IF (SOLID(IC)) CYCLE
            IW = WALL_INDEX(IC,-IOR)
            AREA = AW(IW)
            IF (PRESSURE_ZONE_WALL(IW) /= NODE_ZONE) THEN
               IF (NODE_ZONE == -1) THEN
                  NODE_ZONE = PRESSURE_ZONE_WALL(IW)
               ELSE
                  WRITE(MESSAGE,'(A,I2)') 'ERROR: VENT for a NODE must lie with a single pressure zone. Node: ',NN
                  CALL SHUTDOWN(MESSAGE)
               ENDIF
            ENDIF
            RHO_SUM = RHO_SUM + RHOP(II,JJ,KK)*AW(IW)
            IF (N_SPECIES > 0) THEN
               YY_SUM = YY_SUM+YY(II,JJ,KK,:)*RHO(II,JJ,KK)*AREA
               YY_N = YY(II,JJ,KK,:)
               CALL GET_AVERAGE_SPECIFIC_HEAT(YY_N,H_G,MIN(5000,NINT(TMP(II,JJ,KK))))
               H_G = H_G * TMP(II,JJ,KK)
            ELSE
               H_G = TMP(II,JJ,KK)*Y2CP_C(MIN(5000,NINT(TMP(II,JJ,KK))))
            ENDIF
            H_SUM = H_SUM+H_G*AREA*RHOP(II,JJ,KK)
            TMP_SUM = TMP_SUM + TMP(II,JJ,KK)*AREA
            VEL2  = 0.25_EB*( (UP(II,JJ,KK)+UP(II-1,JJ,KK))**2+(VP(II,JJ,KK)+VP(II,JJ-1,KK))**2+(WP(II,JJ,KK)+WP(II,JJ,KK-1))**2)
            !P_SUM = P_SUM + PBARP(KK,PRESSURE_ZONE(II,JJ,KK))*AREA
            SELECT CASE (IOR)
               CASE (3)
                  P_AVE = 0.5_EB*(PBARP(KK-1,PRESSURE_ZONE(II,JJ,KK-1))+PBARP(KK,PRESSURE_ZONE(II,JJ,KK)))
               CASE (-3)
                  P_AVE = 0.5_EB*(PBARP(KK,PRESSURE_ZONE(II,JJ,KK))+PBARP(KK+1,PRESSURE_ZONE(II,JJ,KK+1)))
               CASE DEFAULT
                  P_AVE = PBARP(KK,PRESSURE_ZONE(II,JJ,KK))
            END SELECT
            P_SUM = P_SUM + (P_AVE+RHOP(II,JJ,KK)*(H(II,JJ,KK)-.5_EB*VEL2))*AREA
            AREA_SUM = AREA_SUM + AREA
         ENDDO
      ENDDO
   ENDDO
   IF (NODE_ZONE /= -1) DN%ZONE_INDEX = NODE_ZONE
   DN%P      = P_SUM / AREA_SUM
   TMP_SUM   = TMP_SUM / AREA_SUM
   IF (N_SPECIES>0) DN%YY_V(:) = YY_SUM / RHO_SUM
   H_SUM = H_SUM / RHO_SUM
   DO
      IF (N_SPECIES > 0) THEN
         YY_N = DN%YY_V(:)
         CALL GET_AVERAGE_SPECIFIC_HEAT(YY_N,CP,MIN(5000,NINT(TMP_SUM)))      
      ELSE
         CP = Y2CP_C(MIN(5000,NINT(TMP_SUM)))
      ENDIF
      TMP_NEW = H_SUM / CP
      IF (ABS(TMP_SUM - TMP_NEW) < ZERO_P .OR. ABS(TMP_SUM - TMP_NEW)/TMP_NEW < 0.0005_EB) EXIT
      TMP_SUM = TMP_NEW
   ENDDO
   DN%CP_V  = CP
   DN%TMP_V = TMP_SUM
   IF (N_SPECIES>0) THEN
      YY_N = DN%YY_V(:)
      CALL GET_SPECIFIC_GAS_CONSTANT(DN%YY_V(:),DN%RSUM_V)
   ENDIF
   DN%RHO_V  = RHO_SUM / AREA_SUM
ENDDO

TUSED(16,:)=TUSED(16,:)+SECOND()-TNOW

END SUBROUTINE HVAC_BC_IN

SUBROUTINE DETERMINE_FIXED_ELEMENTS
INTEGER:: NN,ND, COUNTER,DUCT_INDEX
REAL(EB) :: VOLUME_FLOW,ABS_VOLUME_FLOW
LOGICAL :: CHANGE
TYPE(DUCT_TYPE), POINTER :: DU=>NULL()
TYPE(DUCTNODE_TYPE), POINTER :: DN=>NULL()!,DN2=>NULL()

DUCTNODE%FIXED = .FALSE.
DUCT%FIXED = .FALSE.
DUCT%VOLUME_FLOW = 0._EB
CHANGE = .TRUE.
FIXED_LOOP: DO WHILE (CHANGE)
   CHANGE = .FALSE.
   NODE_LOOP: DO NN = 1, N_DUCTNODES
      DN=>DUCTNODE(NN)
      IF( DN%FIXED) CYCLE NODE_LOOP
      IF (DN%VENT_INDEX > 0 .OR. DN%AMBIENT) THEN
         DN%FIXED = .TRUE.
         CHANGE = .TRUE.
         DU=>DUCT(DN%DUCT_INDEX(1))
         IF (DU%VOLUME_FLOW_INITIAL<1.E6_EB) THEN
            DU%VOLUME_FLOW = DU%VOLUME_FLOW_INITIAL
            IF(DU%AREA > ZERO_P) DU%VEL = DU%VOLUME_FLOW/DU%AREA
            DU%FIXED = .TRUE.
         ENDIF
         CYCLE
      ENDIF
      COUNTER = 0
      DO ND=1,DN%N_DUCTS
         DU=>DUCT(DN%DUCT_INDEX(ND))
         IF (DU%VOLUME_FLOW_INITIAL<1.E6_EB) THEN
            DU%VOLUME_FLOW = DU%VOLUME_FLOW_INITIAL
            DU%FIXED = .TRUE.
         ENDIF
         IF (DU%FIXED .OR. DU%AREA < ZERO_P) THEN
            COUNTER = COUNTER + 1
         ENDIF
      ENDDO
      IF (COUNTER >= DN%N_DUCTS - 1) THEN
         VOLUME_FLOW = 0._EB
         ABS_VOLUME_FLOW = 0._EB      
         DO ND=1,DN%N_DUCTS
            DU=>DUCT(DN%DUCT_INDEX(ND))         
            IF (DU%FIXED) THEN
               IF (DU%AREA > ZERO_P) THEN
                  VOLUME_FLOW = VOLUME_FLOW + DN%DIR(ND) * DU%VOLUME_FLOW
                  ABS_VOLUME_FLOW = ABS_VOLUME_FLOW + ABS(DU%VOLUME_FLOW)
               ENDIF
            ELSE
               DUCT_INDEX = ND
            ENDIF            
         ENDDO
         IF (COUNTER==DN%N_DUCTS) THEN
            IF (ABS(VOLUME_FLOW)/ABS_VOLUME_FLOW > 1.E-6_EB) THEN
               WRITE(MESSAGE,'(A,A)') 'ERROR: Volume flow does not sum to zero for DUCTNODE = ',TRIM(DN%ID)
               CALL SHUTDOWN(MESSAGE)
            ENDIF
         ELSE
            DU=>DUCT(DN%DUCT_INDEX(DUCT_INDEX))
            DU%VOLUME_FLOW = -DN%DIR(DUCT_INDEX) * VOLUME_FLOW
            DU%VEL = DU%VOLUME_FLOW/DU%AREA
            DU%FIXED = .TRUE.
         ENDIF
         DN%FIXED = .TRUE.
         CHANGE = .TRUE.
      ELSEIF(COUNTER==DN%N_DUCTS) THEN
         DN%FIXED = .TRUE.
         CHANGE = .TRUE.         
      ENDIF
   END DO NODE_LOOP
END DO FIXED_LOOP

END SUBROUTINE DETERMINE_FIXED_ELEMENTS


SUBROUTINE FIND_NETWORKS(CHANGEIN)
INTEGER:: NZ,NN,ND,DUCT_COUNTER(N_DUCTS),NODE_COUNTER(N_DUCTNODES),COUNTER,COUNTER2,NETWORK_COUNTER(N_DUCTS),&
          ZONE_COUNTER(N_ZONE)
INTEGER, DIMENSION(:), ALLOCATABLE :: NETWORK_DCOUNTER,NETWORK_NCOUNTER
LOGICAL, INTENT(INOUT) :: CHANGEIN
LOGICAL :: CHANGE!,ZONE_CHANGE
TYPE(DUCT_TYPE), POINTER :: DU=>NULL()
!TYPE(DUCTNODE_TYPE), POINTER :: DN=>NULL(),DN2=>NULL()

CHANGE = CHANGEIN

IF (N_ZONE > 0) ZONE_COUNTER = 0

DO ND = 1, N_DUCTS
   DU => DUCT(ND)
   IF(.NOT. DU%DAMPER) CYCLE
   IF (DU%DEVC_INDEX > 0) THEN
      IF (DEVICE(DU%DEVC_INDEX)%CURRENT_STATE .NEQV. DU%DAMPER_OPEN) THEN
            DU%DAMPER_OPEN = DEVICE(DU%DEVC_INDEX)%CURRENT_STATE
            CHANGE = .TRUE.
            IF (DU%DAMPER_OPEN) THEN
               DU%AREA = DU%AREA_INITIAL
            ELSE
               DU%AREA = 0._EB
               DU%VEL  = 0._EB
            ENDIF
      ENDIF
   ELSE
      IF (CONTROL(DU%CTRL_INDEX)%CURRENT_STATE .NEQV. DU%DAMPER_OPEN) THEN
            DU%DAMPER_OPEN = CONTROL(DU%CTRL_INDEX)%CURRENT_STATE
            CHANGE = .TRUE.
            IF (DU%DAMPER_OPEN) THEN
               DU%AREA = DU%AREA_INITIAL
            ELSE
               DU%AREA = 0._EB
               DU%VEL  = 0._EB
            ENDIF
      ENDIF
   ENDIF
ENDDO
IF (CHANGE) THEN
   IF (ALLOCATED(NETWORK)) DEALLOCATE(NETWORK)
   DO NN = 1, N_DUCTNODES
      NZ = DUCTNODE(NN)%ZONE_INDEX
      IF (NZ>=1) THEN
         NODE_COUNTER(NN) = NZ
         ZONE_COUNTER(NZ) = ZONE_COUNTER(NZ) + 1
      ELSE
         NODE_COUNTER(NN) = NN+N_ZONE
      ENDIF
   ENDDO 
   IF (N_ZONE > 0) THEN
      DO NZ = 1, N_ZONE
       IF (ASSOCIATED(P_ZONE(NZ)%NODE_INDEX)) DEALLOCATE(P_ZONE(NZ)%NODE_INDEX)
         ALLOCATE(P_ZONE(NZ)%NODE_INDEX(ZONE_COUNTER(NZ)))      
         P_ZONE(NZ)%N_DUCTNODES = ZONE_COUNTER(NZ)         
         COUNTER = 1
         DO NN = 1,N_DUCTNODES
            IF (DUCTNODE(NN)%ZONE_INDEX == NZ) THEN
               P_ZONE(NZ)%NODE_INDEX(COUNTER)=NN
               COUNTER = COUNTER + 1
            ENDIF
         ENDDO   
      ENDDO      
   ENDIF  
   CHANGE = .TRUE.
   DO WHILE (CHANGE)
      CHANGE = .FALSE.
      DO ND = 1, N_DUCTS
         DU => DUCT(ND)
         IF (NODE_COUNTER(DU%NODE_INDEX(1)) /= NODE_COUNTER(DU%NODE_INDEX(2))) THEN
            CHANGE = .TRUE.
            COUNTER = MIN(NODE_COUNTER(DU%NODE_INDEX(1)),NODE_COUNTER(DU%NODE_INDEX(2)))
            DUCT_COUNTER(ND) = COUNTER
            NODE_COUNTER(DU%NODE_INDEX(1)) = COUNTER
            NODE_COUNTER(DU%NODE_INDEX(2)) = COUNTER
         ELSE
            DUCT_COUNTER(ND) = NODE_COUNTER(DU%NODE_INDEX(1))
         ENDIF
      ENDDO 
      IF (N_ZONE > 0) THEN   
         DO NZ = 1, N_ZONE
            COUNTER = 1
            COUNTER2 = 1
            DO NN = 1, P_ZONE(NZ)%N_DUCTNODES
               IF (NN==1) THEN
                  COUNTER = NODE_COUNTER(P_ZONE(NZ)%NODE_INDEX(NN))
                  COUNTER2 = COUNTER
               ELSE
                  IF (COUNTER /= NODE_COUNTER(P_ZONE(NZ)%NODE_INDEX(NN))) &
                     COUNTER2 = MAX(COUNTER2,NODE_COUNTER(P_ZONE(NZ)%NODE_INDEX(NN)))                     
                     COUNTER  = MIN(COUNTER,NODE_COUNTER(P_ZONE(NZ)%NODE_INDEX(NN)))                     
               ENDIF               
            ENDDO
            IF (COUNTER /= COUNTER2) THEN
               CHANGE = .TRUE.
               DO NN = 1, P_ZONE(NZ)%N_DUCTNODES   
                  NODE_COUNTER(P_ZONE(NZ)%NODE_INDEX(NN)) = COUNTER
               ENDDO
            ENDIF
         ENDDO
      ENDIF
   END DO
   NETWORK_COUNTER = 0
   DO ND = 1, N_DUCTS
      NETWORK_COUNTER(DUCT_COUNTER(ND))=1
   ENDDO

   N_NETWORKS = SUM(NETWORK_COUNTER)
   ALLOCATE(NETWORK(N_NETWORKS))
   NETWORK%N_DUCTS=0
   NETWORK%N_DUCTNODES=0
   ALLOCATE(NETWORK_DCOUNTER(N_NETWORKS))
   NETWORK_DCOUNTER=0
   ALLOCATE(NETWORK_NCOUNTER(N_NETWORKS))
   NETWORK_NCOUNTER=0
   COUNTER = 0
   DO ND = 1, N_DUCTS
      IF(NETWORK_COUNTER(ND)==1) THEN
         COUNTER = COUNTER + 1
         NETWORK_COUNTER(ND) = COUNTER
      ENDIF
   ENDDO   
   DO ND = 1, N_DUCTS
      NETWORK(NETWORK_COUNTER(DUCT_COUNTER(ND)))%N_DUCTS = NETWORK(NETWORK_COUNTER(DUCT_COUNTER(ND)))%N_DUCTS + 1
   ENDDO
   DO NN = 1, N_DUCTNODES
      NETWORK(NETWORK_COUNTER(NODE_COUNTER(NN)))%N_DUCTNODES = NETWORK(NETWORK_COUNTER(NODE_COUNTER(NN)))%N_DUCTNODES + 1
   ENDDO
   DO NN = 1, N_NETWORKS
      ALLOCATE(NETWORK(NN)%DUCT_INDEX(NETWORK(NN)%N_DUCTS))
      ALLOCATE(NETWORK(NN)%NODE_INDEX(NETWORK(NN)%N_DUCTNODES))
      ALLOCATE(NETWORK(NN)%MATRIX_INDEX(NETWORK(NN)%N_DUCTS+NETWORK(NN)%N_DUCTNODES))
      NETWORK(NN)%MATRIX_INDEX = 0
   ENDDO
   DO ND = 1, N_DUCTS
      NETWORK_DCOUNTER(NETWORK_COUNTER(DUCT_COUNTER(ND))) = NETWORK_DCOUNTER(NETWORK_COUNTER(DUCT_COUNTER(ND))) + 1
      NETWORK(NETWORK_COUNTER(DUCT_COUNTER(ND)))%DUCT_INDEX(NETWORK_DCOUNTER(NETWORK_COUNTER(DUCT_COUNTER(ND)))) = ND
      DUCT_NE(ND) = NETWORK_DCOUNTER(NETWORK_COUNTER(DUCT_COUNTER(ND)))
   ENDDO
   DO NN = 1, N_DUCTNODES
      NETWORK_NCOUNTER(NETWORK_COUNTER(NODE_COUNTER(NN))) = NETWORK_NCOUNTER(NETWORK_COUNTER(NODE_COUNTER(NN))) + 1
      NETWORK(NETWORK_COUNTER(NODE_COUNTER(NN)))%NODE_INDEX(NETWORK_NCOUNTER(NETWORK_COUNTER(NODE_COUNTER(NN)))) = NN
      DUCTNODE_NE(NN) = NETWORK_NCOUNTER(NETWORK_COUNTER(NODE_COUNTER(NN)))
   ENDDO
   DEALLOCATE(NETWORK_DCOUNTER)
   DEALLOCATE(NETWORK_NCOUNTER)
   CALL DETERMINE_FIXED_ELEMENTS
   CALL SETUP_SOLUTION_POINTERS
   !WRITE(*,*) N_NETWORKS,NETWORK(1)%N_DUCTS,NETWORK(1)%N_DUCTNODES
   !WRITE(*,*) NETWORK(1)%DUCT_INDEX
   !WRITE(*,*) NETWORK(1)%NODE_INDEX  
   !WRITE(*,*) NETWORK(1)%N_MATRIX 
   !WRITE(*,*) NETWORK(1)%MATRIX_INDEX
   !WRITE(*,*) DUCT_NE
   !WRITE(*,*) DUCTNODE_NE
ENDIF

END SUBROUTINE FIND_NETWORKS

SUBROUTINE SETUP_SOLUTION_POINTERS
INTEGER:: NNE,NN,ND,COUNTER
TYPE(DUCT_TYPE), POINTER :: DU=>NULL()
TYPE(DUCTNODE_TYPE), POINTER :: DN=>NULL()
TYPE(NETWORK_TYPE), POINTER :: NE=>NULL()

DO NNE = 1,N_NETWORKS
   COUNTER = 0
   NE => NETWORK(NNE)   
   DO ND=1,NE%N_DUCTS
      DU=>DUCT(NE%DUCT_INDEX(ND))
      IF (DU%FIXED .OR. DU%AREA<ZERO_P) CYCLE
      COUNTER = COUNTER + 1
      NE%MATRIX_INDEX(ND)=COUNTER
   ENDDO
   DO NN=1,NE%N_DUCTNODES
      DN=>DUCTNODE(NE%NODE_INDEX(NN))      
      IF (DN%FIXED .OR. DN%VENT_INDEX>0) CYCLE
      COUNTER = COUNTER + 1
      NE%MATRIX_INDEX(NE%N_DUCTS+NN)=COUNTER
   ENDDO
   NE%N_MATRIX=COUNTER
ENDDO

END SUBROUTINE SETUP_SOLUTION_POINTERS


SUBROUTINE UPDATE_LOSS(NNE)
REAL(EB) :: FRICTION_FACTOR,LOSS_SUM
INTEGER, INTENT(IN) :: NNE
INTEGER :: ND,ND2, NN
TYPE(DUCT_TYPE), POINTER :: DU=>NULL()
TYPE(DUCTNODE_TYPE), POINTER :: DN=>NULL()
TYPE(NETWORK_TYPE), POINTER :: NE=>NULL()

NE => NETWORK(NNE)
!!IF(N_FILTERS > 0) CALL FILTER_UPDATE(NNE)
DO ND = 1, NE%N_DUCTS
   DUCT(NE%DUCT_INDEX(ND))%TOTAL_LOSS=0._EB
ENDDO
NODELOOP : DO NN=1,NE%N_DUCTNODES
  DN => DUCTNODE(NE%NODE_INDEX(NN))
  DN%BRANCH_LOSS = 0._EB
  NODECLASS: IF (DN%FILTER_INDEX > 0) THEN
     IF(DUCT(DN%DUCT_INDEX(1))%VEL(GUESS)*DN%DIR(1) < 0._EB) THEN
        DN%BRANCH_LOSS(2) = DN%FILTER_LOSS
     ELSE
        DN%BRANCH_LOSS(1) = DN%FILTER_LOSS     
     ENDIF
  ELSEIF(DN%VENT_INDEX > 0) THEN NODECLASS
     IF(DUCT(DN%DUCT_INDEX(1))%FIXED) CYCLE
     IF(DUCT(DN%DUCT_INDEX(1))%VEL(GUESS)*DN%DIR(1) > 0._EB) THEN
        DN%BRANCH_LOSS(1) = DN%LOSS_ARRAY(1,2)
     ELSE
        DN%BRANCH_LOSS(1) = DN%LOSS_ARRAY(2,1)
     ENDIF     
  ELSE NODECLASS
     DO ND=1,DN%N_DUCTS
        DU => DUCT(DN%DUCT_INDEX(ND))
        IF (DU%VEL(GUESS)*DN%DIR(ND) > 0._EB) THEN
           DN%BRANCH_LOSS(ND) = ABS(DU%VEL(GUESS))*DU%AREA*DU%RHO_D
        ELSE
           DN%BRANCH_LOSS(ND) = 0._EB
        ENDIF
     ENDDO
     IF (SUM(DN%BRANCH_LOSS) < ZERO_P) THEN
        DN%BRANCH_LOSS = 0.5_EB * SUM(DN%LOSS_ARRAY)/DN%N_DUCTS
     ELSE
        DN%BRANCH_LOSS = DN%BRANCH_LOSS / SUM(DN%BRANCH_LOSS)
        DO ND=1,DN%N_DUCTS
           IF (DN%BRANCH_LOSS(ND)> ZERO_P) CYCLE
           LOSS_SUM = 0._EB
           DO ND2=1,DN%N_DUCTS
             IF(DN%BRANCH_LOSS(ND2) < ZERO_P) CYCLE
             LOSS_SUM = LOSS_SUM + DN%BRANCH_LOSS(ND2)*DN%LOSS_ARRAY(ND,ND2)
           ENDDO
           DUCT(DN%DUCT_INDEX(ND))%TOTAL_LOSS = LOSS_SUM
        ENDDO
     ENDIF
  ENDIF NODECLASS
ENDDO NODELOOP
DO ND = 1, NE%N_DUCTS
   DU => DUCT(NE%DUCT_INDEX(ND))
   IF (DU%ROUGHNESS > ZERO_P) THEN
      FRICTION_FACTOR = COMPUTE_FRICTION_FACTOR(DU%RHO_D,DU%TMP_D,ABS(DU%VEL(GUESS)),DU%DIAMETER,DU%ROUGHNESS)
   ELSE
      FRICTION_FACTOR = 0._EB
   ENDIF
   IF (DU%VEL(GUESS)>0._EB) THEN
      DU%TOTAL_LOSS = DU%LOSS(1)+DU%TOTAL_LOSS
   ELSEIF (DU%VEL(GUESS)<0._EB) THEN
      DU%TOTAL_LOSS = DU%LOSS(2)+DU%TOTAL_LOSS
   ELSE
      DU%TOTAL_LOSS = DU%TOTAL_LOSS+0.5_EB*(DU%LOSS(1)+DU%LOSS(2))
   ENDIF
   DU%TOTAL_LOSS = DU%TOTAL_LOSS + DU%LENGTH/DU%DIAMETER*FRICTION_FACTOR
   IF (DU%FAN_INDEX>0) THEN
      IF(.NOT. DU%FAN_OPERATING) DU%TOTAL_LOSS = DU%TOTAL_LOSS + FAN(DU%FAN_INDEX)%OFF_LOSS
   ENDIF
ENDDO

END SUBROUTINE UPDATE_LOSS


REAL(EB) FUNCTION COMPUTE_FRICTION_FACTOR(RHO,TEMPIN,VEL,DIAM,ROUGHNESS)
REAL(EB), INTENT(IN) :: RHO,TEMPIN,VEL,DIAM,ROUGHNESS
REAL(EB) :: EOD,RE_D,VISCOSITY,TEMP

TEMP = MAX(MIN(TEMPIN,2400._EB),200._EB)
!Visocisty power interpolation of 200 K - 2400 K, 200 K step,  J.P. Holman, Heat Transfer, 7th Ed
VISCOSITY = 3.9424E-7_EB*TEMP**0.67311_EB
RE_D = MAX(100._EB,RHO*DIAM*VEL/VISCOSITY)
EOD = ROUGHNESS / DIAM
COMPUTE_FRICTION_FACTOR = LOG10(6.9_EB/RE_D+(EOD/3.7_EB)**1.11_EB)
COMPUTE_FRICTION_FACTOR = LOG10(EOD/3.7_EB-4.518_EB/RE_D*COMPUTE_FRICTION_FACTOR)
COMPUTE_FRICTION_FACTOR = (-0.5_EB/COMPUTE_FRICTION_FACTOR)**2

RETURN

END FUNCTION COMPUTE_FRICTION_FACTOR


!SUBROUTINE FILTER_UPDATE(NNE)
!INTEGER, INTENT(IN) :: NNE
!
!RETURN
!
!END SUBROUTINE FILTER_UPDATE


SUBROUTINE SET_GUESS(NNE,ITER)
INTEGER, INTENT(IN) :: NNE,ITER
INTEGER :: ND
REAL(EB) :: VEL_TMP
TYPE(DUCT_TYPE),POINTER :: DU=>NULL()
TYPE(NETWORK_TYPE), POINTER :: NE=>NULL()

NE => NETWORK(NNE)

DO ND = 1,NE%N_DUCTS
     DU => DUCT(NE%DUCT_INDEX(ND))
      IF (ITER == 0) THEN
         DU%VEL(GUESS)     = DU%VEL(OLD)
         DU%VEL(NEW)       = DU%VEL(OLD)
         DU%VEL(PREVIOUS)  = DU%VEL(OLD)  
      ELSE
         VEL_TMP = DU%VEL(NEW)
         IF (ABS(DU%VEL(GUESS)) < ZERO_P) THEN
            DU%VEL(GUESS)     = DU%VEL(NEW)         
            DU%VEL(NEW)       = DU%VEL(PREVIOUS)
         ELSEIF (SIGN(1._EB,DU%VEL(PREVIOUS))==SIGN(1._EB,DU%VEL(NEW))) THEN
            DU%VEL(GUESS)     = DU%VEL(NEW)         
            DU%VEL(NEW)       = 0.6_EB*DU%VEL(PREVIOUS)+0.4_EB*DU%VEL(GUESS)
         ELSE
            DU%VEL(GUESS)     = DU%VEL(NEW)         
            DU%VEL(NEW)       = 0._EB
         ENDIF
         DU%VEL(PREVIOUS) = VEL_TMP
      ENDIF
   ENDDO

END SUBROUTINE SET_GUESS


SUBROUTINE SET_DONOR(NNE,T)
USE MATH_FUNCTIONS, ONLY : EVALUATE_RAMP
INTEGER :: ND
REAL(EB), INTENT(IN) :: T
INTEGER, INTENT(IN) :: NNE
REAL(EB) :: RHOLAST,TMPLAST,FVAL,OMFVAL,FLOW_FAC,ITERFRAC
TYPE(DUCT_TYPE), POINTER :: DU=>NULL()
TYPE(DUCTNODE_TYPE), POINTER :: DN=>NULL()
TYPE(NETWORK_TYPE), POINTER :: NE=>NULL()

NE => NETWORK(NNE)
ITERFRAC = REAL(ITER,EB)/REAL(ITER_MAX,EB)
FVAL = MIN(1._EB,MAX(0._EB,(ITERFRAC-ONTH))/ONTH,1._EB)
OMFVAL = 1._EB - FVAL

DUCTLOOP: DO ND=1,NE%N_DUCTS
   DU=>DUCT(NE%DUCT_INDEX(ND))
   IF (DU%AREA < ZERO_P) CYCLE DUCTLOOP
   RHOLAST = DU%RHO_D   
   TMPLAST = DU%TMP_D
   IF (DU%VEL(PREVIOUS)>ZERO_P) THEN
      DN=>DUCTNODE(DU%NODE_INDEX(1))
      IF(DN%VENT_INDEX > 0 .OR. DN%AMBIENT) THEN
         DN%RHO= DN%RHO_V
         DN%TMP= DN%TMP_V
         IF (N_SPECIES>0) DN%YY = DN%YY_V
         DN%RSUM = DN%RSUM_V
         DN%CP = DN%CP_V                  
      ENDIF
   ELSE
      DN=>DUCTNODE(DU%NODE_INDEX(2))     
      IF(DN%VENT_INDEX > 0 .OR. DN%AMBIENT) THEN
         DN%RHO= DN%RHO_V
         DN%TMP= DN%TMP_V
         IF (N_SPECIES>0) DN%YY = DN%YY_V
         DN%RSUM = DN%RSUM_V         
         DN%CP = DN%CP_V                  
      ENDIF
   ENDIF
   DU%RHO_D = DN%RHO
   DU%TMP_D = DN%TMP
   DU%RSUM_D = DN%RSUM
   IF (ITERFRAC > ONTH) THEN
      DU%RHO_D = FVAL*RHOLAST + OMFVAL*DU%RHO_D
      DU%TMP_D = FVAL*TMPLAST + OMFVAL*DU%TMP_D
   ENDIF   
   IF (DU%FIXED) THEN
      IF (DU%RAMP_INDEX > 0) THEN
         FLOW_FAC = EVALUATE_RAMP(T,1._EB,DU%RAMP_INDEX)
         DU%VOLUME_FLOW = DU%VOLUME_FLOW_INITIAL * FLOW_FAC
      ENDIF
      DU%VEL(GUESS) = DU%VOLUME_FLOW/DU%AREA
      DU%VEL(PREVIOUS) = DU%VEL(GUESS)
      DU%VEL(OLD)   = DU%VEL(GUESS)
      DU%VEL(NEW)   = DU%VEL(GUESS)
   ENDIF   
ENDDO DUCTLOOP

END SUBROUTINE SET_DONOR


SUBROUTINE CONVERGENCE_CHECK(ITER,NNE)
INTEGER, INTENT(INOUT) :: ITER
INTEGER, INTENT(IN) :: NNE
INTEGER :: NN, ND
LOGICAL :: CONVERGED
REAL(EB) :: MSUM,MTOT,MFLOW,VEL
TYPE(NETWORK_TYPE), POINTER :: NE=>NULL()
TYPE(DUCT_TYPE), POINTER :: DU=>NULL()
TYPE(DUCTNODE_TYPE), POINTER :: DN=>NULL()

NE => NETWORK(NNE)
CONVERGED = .TRUE.
!Check duct velocity convergence
DO ND=1,NE%N_DUCTS

   DU => DUCT(NE%DUCT_INDEX(ND))
   IF (DU%FIXED .OR. DU%AREA < ZERO_P) CYCLE
   IF (ABS(DU%VEL(PREVIOUS)) < 1.E-5_EB .AND. ABS(DU%VEL(NEW)) < 1.E-5_EB) CYCLE
   IF (DU%VEL(PREVIOUS) < 0._EB .EQV. DU%VEL(NEW) < 0._EB) THEN
        IF (ABS(DU%VEL(PREVIOUS))<ZERO_P) THEN
           CONVERGED = .FALSE.
           EXIT
        ELSE
           IF (ABS(1._EB-DU%VEL(NEW)/DU%VEL(PREVIOUS)) > 0.05_EB) CYCLE
        ENDIF
   ENDIF
   CONVERGED = .FALSE.
ENDDO

IF (.NOT. CONVERGED) RETURN

!Check node mass conservation convergence
DO NN=1,NE%N_DUCTNODES
   DN => DUCTNODE(NE%NODE_INDEX(NN))
   IF (DN%FIXED) CYCLE
   MSUM = 0._EB
   MTOT = 0._EB
   DO ND=1,DN%N_DUCTS
      DU=>DUCT(DN%DUCT_INDEX(ND))
      IF (ABS(DU%VEL(NEW))<1.E-7_EB) THEN
         VEL = 0._EB
      ELSE
         VEL = DU%VEL(NEW)
      ENDIF
      MFLOW = DN%DIR(ND)*VEL*DU%RHO_D*DU%AREA
      MSUM = MSUM + MFLOW
      MTOT = MTOT + ABS(MFLOW)
   ENDDO
   IF(ABS(MSUM)< 1.E-6 * MTOT .OR. MTOT < ZERO_P) CYCLE
   CONVERGED = .FALSE.
ENDDO

IF (CONVERGED) ITER=ITER_MAX

END SUBROUTINE CONVERGENCE_CHECK


SUBROUTINE GET_REV_HVAC(MODULE_REV,MODULE_DATE)
INTEGER,INTENT(INOUT) :: MODULE_REV
CHARACTER(255),INTENT(INOUT) :: MODULE_DATE
WRITE(MODULE_DATE,'(A)') hvacrev(INDEX(hvacrev,':')+1:LEN_TRIM(hvacrev)-2)
READ (MODULE_DATE,'(I5)') MODULE_REV
WRITE(MODULE_DATE,'(A)') hvacdate
END SUBROUTINE GET_REV_HVAC
 
END MODULE HVAC_ROUTINES
